arXiv:1504.00214v2 [physics.plasm-ph] 2 Apr 2015 


Adjoint Monte Carlo calculation of charged plasma particle flux to wall 


Simppa Akaslompolo 

Aalto University, Finland 


Abstract 

This manuscript describes an adjoint/reverse Monte Carlo method to calculate the flux of charged plasma particles to 
the wall of e.g. a tokamak. Two applications are described: a fusion product activation probe and a neutral beam 
injection prompt loss measurement with a fast ion loss diagnostic. In both cases, the collisions of the particles with the 
background plasma can be omitted. 


1. Introduction 

In the research of magnetically confined fusion, there is 
a need to study the behaviour of minority particles in the 
plasma. A typical example is the calculation of the density 
distribution of slowing down fusion alpha particles in a 
tokamak reactor. This characterises the ions confined in 
the plasma. Another interesting feature are the ions that 
escape the plasma and hit the walls. The escaping ions 
can damage the wall, but on the other hand, the plasma 
can be studied by measuring the escaping ion flux. 

The usual tool for this kind of studies is a Monte Carlo 
particle following code that draws an ensemble of mark¬ 
ers from a source distribution, such as fusion reactivity. 
The trajectory of each marker is then followed accord¬ 
ing to the equations of motion derived from the Lorentz 
force, while applying a stochastic diffusion operator to 
account for collisions of the marker with the background 
plasma. The author is a member of the group using and 
developing the ASCOT code fT], so it is used as the model 
case. Other similar tools include at least OFMC [J], SPI¬ 
RAL Ell SPOT a and LOCUST-GPU HH. These codes 
work well for the calculation of the distribution function, 
and also for the escaping flux to the large features of the 
wall. However, if we are interested in the flux to a small 
detail of the wall, the target, the calculation quickly be¬ 
comes highly inefficient: from the large initial ensemble 
only a vanishing fraction hits the target. 

The adjoint Monte Carlo integration is based on the 
idea of swapping the roles of the relatively large sources 
and tiny targets in the calculation. 
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2. Proposed adjoint calculation method 

In our case the source S is a source of particles within 
the plasma, such as fusion reactivity in activation probe 
simulations |[6|. The target could be any small surface of 
the wall where the flux of particles is of interest. Since 
the target is small and the source is large, very few mark¬ 
ers will find their way to the target, so most of them are 
“wasted”. In the adjoint method, however, markers start 
backward in time from the target and are much more likely 
to reach the source. 

The idea of the method is to check how much flux 
(how many particles per second) arrives at the target through 
all possible final states of the particles that have reached 
the target. Markers are then followed backwards in time 
from the final states and the sum (integral) over all the 
source that would have contributed to the final states is 
gathered. The trick is to only follow particles backwards 
in time until they have experienced their last collision be¬ 
fore hitting the target. The method could be called final 
collision integral. 

All possible final states are included by forming a Mon¬ 
te Carlo integral of them. This is the basis for the initial¬ 
isation. In other words, N markers are launched from the 
wall element so that each one carries the weight equiva¬ 
lent of (Aip • Av • Cvr) /A, where Ayy is the surface area of 
the target, Av = v^ax - is the velocity range of the 
final state and dCvv is the solid angle of the possible final 
velocity directions. For a planar target, the final veloc¬ 
ity distribution is isotropic in the half-space open to the 
plasma (Ow = In), and the final locations are uniformly 
distributed on the surface of the target. (Additionally, im¬ 
portance sampling could be applied, but for the sake of 
simplicity identical weights are used here.) 
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In the time-reversed calculation, the full gyromotion 
of the particles is followed without collisions. At each 
time step, a certain phase space volume’s worth {dV ■ dO • 
dv) of source S is added to the flux integral. Naturally, 
particles on the same trajectory with exactly the veloc¬ 
ity of the time-reversed marker would hit the target, but 
particles with nearly same location and nearly the same 
velocity would reach the target as well. The above vague 
statement is next refined for use in calculating the phase 
space volume. 

Figure [T] illustrates how the phase space volume ele¬ 
ment dV ■ dfl • dv are related to the Monte Carlo integral 
weights: 

. . Aw • Av • flw 

dAw ■ dv • dOiy =-—- (1) 


The key idea behind the concept is considering that the 
target wall surface element dAw is visible in the solid an¬ 
gle dfl at the particle location. On the other hand, the solid 
angle element dOvv and the (reverse) time step dt define a 
cylinder-like volume dV. The calculation is based on the 
general solid angle formula Q. = Ajp-. The distance r = vt 
can be calculated along the particle trajectory instead of a 
straight line, as Liouville’s theorem suggests that a curved 
path along the particle trajectory will conserve the phase 
space element. We arrive at the following equations: 


dO 

dV 


dAw 

vdt • dQvv • 


( 2 ) 

( 3 ) 


where the = (vt)- part actually cancels out when the 
phase space volume element components dv, dO and dV 
are multiplied. 

The modelling of the collisions of fast particles back¬ 
ward in time is at least non-trivial, if not impossible. This 
is suggested by the forward collisions producing a Max¬ 
wellian distribution regardless of the initial distribution. 
Such collisions increase the entropy of the system or in 
other words remove information. Any reverse calculation 
would need to supply information to the system, which 
seems impossible. 

This problem is averted either by simply ignoring the 
collisions, if possible, or following only collisionless or¬ 
bits backwards in time. The missing collisions are com¬ 
pensated for by including all the collisions in a forward 
Monte Carlo calculation that calculates the source S. The 
adjoint integration phase includes only the path the parti¬ 
cles take after their final collision before hitting the target. 
This is achieved by weighting each contribution to the in¬ 
tegral with the probability C of the marker travelling to 
the target from the source location R without a collision. 




Figure 1: Basis for reverse in time integration 
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The integral of particle flux to the wall elements is 
then: 


F - 


^ fill 
ffU 

AiyAvilvT 


S(R,v)C(R,v)dFdndv (4) 


S (R, v)C(R, v)vdtdQH,dvdAiy(5) 


where v depicts the velocity vector of the time-reversed 
marker and is the life time of the marker. 


3. Application for fusion product simulation 

A fusion product activation probe f7]| is a diagnostic probe 
that can measure the number of high energy (~MeV) fu¬ 
sion products impacting on it. It is inserted near the plasma 
in a tokamak. Fusion products activate the samples in the 
probe, and the activity is measured post mortem. 

For collisionless fusion products, the method resem¬ 
bles the calculation of detection efflcience e in [8], but has 
fewer analytic approximations. The method presented in 
this manuscript was indepently developed. 

When simulating the fusion product activation probe 
in ASDEX Upgrade tokamak, we have made certain ap¬ 
proximations: (1) We neglect the collisions of fusion prod¬ 
ucts after birth; in other words C = 1. (2) We use an 
isotropic monoenergetic fusion product source. In reality 
the particles have an anisotropic and relatively wide birth 
energy distribution. The simplified birth distribution al¬ 
lows us to: (a) drop the integral over Av in equation Q 
and (b) consider the solid angle dependency simply by di¬ 
viding by An. This will result in the following source term 
formulation: 

(crv>nin2 

s ^ (6) 

where {crv)n\n 2 is the fusion reactivity in units l/(m^s). 
(In the general case S would be a function of location 
and velocity. Thus it would need to have the units [S ] = 
1 / (s • m/s • m^ • srad^.) 

This will result in the following integral for the parti¬ 
cle flux to the target: 

F- JJ J (7) 

precalculated 

The source term can be separately precalculated into a 
(R, z) grid. The rest of the integral can also be stored as a 
function of R and z, so it too is precalculated and stored 


in a grid, as an adjoint density. This makes it practical to 
quickly change the calculated reactivities or adjoint densi¬ 
ties while keeping the other quantity intact: the final flux 
is calculated as a sum over element-wise multiplied reac¬ 
tivity and time-reversed density. 

The first tests of the above scheme against experimen¬ 
tal results have been reasonably successful @. 

4. Outline of application to NBI prompt losses to a fast 
ion loss detector or Faraday cups 

The the Faraday cups |'9] and the fast ion loss detector 
(FILD) ifTOl are other escaping fast ion diagnostics. The 
calculation method described for an activation probe would 
certainly be adaptable for the FILD, and probably also for 
the Faraday cups, but since the author has no experience 
working with Faraday cups, this discussion focusses on 
the FILD. 

The LILD probe is usually sensitive to neutral beam 
injected ion prompt losses. This means that the fast deu- 
terons (~100keV) hit the walls (or the LILD probe) very 
quickly after being ionized by the plasma. They have little 
to no time to collide with the background plasma. It would 
probably be acceptable to study these losses even neglect¬ 
ing collisions. Lorward Monte Carlo studies have been 
made IfTll . but the simulations require huge resources and 
still cannot provide good statistics in the probe or repro¬ 
duce the detailed probe geometry. Using the adjoint Monte 
Carlo method could provide an efficient calculation of the 
absolutely calibrated flux of prompt ions to the probe, 
possibly with detailed geometry. Time-reversed calcula¬ 
tions are already performed routinely with the GOUR- 
DON code ifT^ to study the final orbits of the measured 
ions, but this is not equivalent to an adjoint calculation. 

When adapting the calculation, the only significant 
modification from the method described for the activation 
probe is changing the source term S from Q to describe 
the NBI ions. Here we present two possible schemes. 

At the moment the initial ensemble of markers for the 
forward Monte Carlo calculation in ASCOT is usually 
created by calculating the ionization locations with the 
Monte Carlo method |[T3l . The straightforward method 
would be to bin the initial ion locations to a multidimen¬ 
sional histogram describing the phase space as a function 
of, e.g., a convenient subspace of {R,(p,z,VR,v^,v^) (ma¬ 
jor radius, toroidal angle and vertical coordinate, and the 
respective velocities) and then use this histogram as S. 

Another possibility would be to utilize the geometrical 
description of the NBI injector model. This is illustrated 
in figure The ionization probability at a given phase 
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space location could possibly be calculated along the way. 
The beam consists of many beamlets which have a source 
point, a direction and typically a (bi)gaussian divergence, 
which means that the beam weakens as a function of an¬ 
gle a shown the figure. The NBI power directed at the 
given phase space location could then be calculated di¬ 
rectly from the NBI model data. However, the beamlet is 
attenuated in the plasma due to more and more atoms get¬ 
ting ionized. This attenuation could conceivably be tabu¬ 
lated prior to the actual adjoint calculation for quick look¬ 
up as a function of distance R from the beamlet source. 
The solid angle df^NB seen by the beamlet source of the 
volume element dV could be calculated as 

lvt\^ 

d^NB - dfliy 1^—j . (8) 

The ionization probability would be ionizations per solid 
angle visible from the injector dQNB> length along beam- 
lets vdt and per second. Finally, the effective solid an¬ 
gle used in the integration should be the “intersection” of 
df^NB and dQ. In other words, the angle /? between the ve¬ 
locities of the neutral and the time-reversed marker must 
be sufficiently small; otherwise, the particle would not hit 
dA]Y. 

The first attempt at prompt loss calculation should prob¬ 
ably be made with the binning method, as the latter method 
seems unnecessarily complicated. 

5. Summary and discussion 

We have described a method to perform adjoint Monte 
Carlo calculation of fusion products. This method has 
been already tested with promising results |j6|. We have 
also outlined a method for neutral beam injection (NBI) 


prompt losses, but this has not been implemented or tested 
yet. 

Future work, would be to test the NBI prompt loss method 
of section The really big deal, however, would be to 
come up with a way to do collisional adjoint calculations. 
This would require either time-reversed collision opera¬ 
tors, or a method to calculate some sort of source term 
for collided particle flux. The section has been written 
with such a source term in mind. Conceivably, the flux of 
particles within the plasma could be calculated with a for¬ 
ward Monte Carlo calculation that includes the collisions. 
This flux could then be “collided” against a suitable differ¬ 
ential (Coulomb) collision/scattering cross section to pro¬ 
duce the sought-after source term. The details of such a 
method remain future work, possibly to be published in a 
future version of this arXiv.org manuscript. 
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